#name : power for 2+ binary/prop
# key : bi.prop.power
# contributor: Shuguang Sun
# --
${1:biprop.power} <- function(p0, p1, p2, n0, n1, n2, alpha = 0.05) {

  z <- qnorm(1 - alpha/2)
  q0 <- 1 - p0
  q1 <- 1 - p1
  q2 <- 1 - p2
  pm1 <- (n0 * p0 + n1 * p1)/(n0 + n1)
  ds1 <- z * sqrt((1/n0 + 1/n1) * pm1 * (1 - pm1))

  pm2 <- (n0 * p0 + n2 * p2)/(n0 + n2)
  ds2 <- z * sqrt((1/n0 + 1/n2) * pm2 * (1 - pm2))

  mean <- c(p1 - p0, p2 - p0)
  sigma <- matrix(c(p0 * q0/n0 + p1 * q1/n1,
                    rep(p0 * q0/n0, 2),
                    p0 * q0/n0 + p2 * q2/n2), nrow = 2)

  pw <- pmvnorm(lower = c(ds1, ds2), upper = Inf,
                mean = mean,
                sigma = sigma, keepAttr = FALSE)

  pw

}

$1($0)